#ifndef GEO_GRID_H
#define GEO_GRID_H

#include <algorithm>
#include <cmath>
#include <vector>

#include "geogram/mesh/mesh.h"
#include "geogram/delaunay/periodic_delaunay_3d.h"
#include "geogram/delaunay/delaunay_2d.h"
#include "geogram/basic/command_line.h"
#include "geogram/basic/command_line_args.h"
#include "geogram/mesh/mesh.h"
#include "geogram/mesh/mesh_geometry.h"
#include "geogram/mesh/mesh_io.h"
#include "geogram/mesh/mesh_remesh.h"

#include "base.h"
#include "utils.h"
#include "geo2d.h"
#include "GRDECLData.h"
#include "kdtree.h"

// using namespace GEO;

// points int a grid layer
class LayerPoints {
public:
    int size[2];
    // points in the layer, points[0] : upper layer, points[1] : lower layer
    std::vector<Point3d> points[2];
};

class GeoDelaunay {
public:
    GEO::PeriodicDelaunay3d* d;
    // GEO::ParallelDelaunay3d* d;
    std::vector<double> flat_points;
    GEO::PeriodicDelaunay3d::IncidentTetrahedra W_;
    bool periodic;

    void SaveToVTK(std::string filename); 
    void get_cell(GEO::index_t v, GEO::ConvexCell& C);
    void Lloyd_iteration();
};

class GeoDelaunay2d {
public:
    struct Point2d {
    };

    GEO::Delaunay2d* d;
    std::vector<double> flat_points;
    int n_points;
    geo2d::Point2d getPointByIndex(int index) {
        return geo2d::Point2d(flat_points[index * 2], flat_points[index * 2 + 1]);
    }
    void SaveToVTK(std::string filename); 
};

class LayerSurface {
public:
    std::vector<Point3d> pts;
    struct Triangle {
        INDEX v[3];
        Point3d p[3];
        Point3d min_brd, max_brd;
        Point2d coord_2d[3];

        Triangle(INDEX v0, INDEX v1, INDEX v2) {
            v[0] = v0;
            v[1] = v1;
            v[2] = v2;
            min_brd = Point3d(1e9, 1e9, 1e9);
            max_brd = Point3d(-1e9, -1e9, -1e9);
        }
        Point3d GetMinBrd() {
            return min_brd;
        }
        Point3d GetMaxBrd() {
            return max_brd;
        }
        Point3d GetCenter() {
            return Point3d((min_brd.x + max_brd.x) / 2.0, (min_brd.y + max_brd.y) / 2.0, (min_brd.z + max_brd.z) / 2.0);
        }
    };
    std::vector<Triangle> tris;
    KDTree<Triangle> kdtree;
    void BuildKDTree() {
        kdtree.BuildTree(tris.data(), tris.size());
    }
    void SaveToVTK(std::string filename);
    // 更新三角形的边界和 2d 坐标
    void UpdateTriRange();
    // 三角形所在平面与 p 所在竖直线的交点 z 坐标
    double TriPlaneIntersection(Triangle tri, Point3d p);
    // 从 p 点发射向上的射线，是否与表面相交 
    int checkInsertWithUpLine(Point3d p);
    // 从 p 点发射向下的射线，是否与表面相交 
    int checkInsertWithDownLine(Point3d p);
};

// 网格延拓信息
class ExpandInfo {
public:
    double expand_top;
    double expand_bottom;
    double expand_X_ratio;
    double expand_Y_ratio;

    Point3d inner_rmin, inner_rmax;
    Point3d outer_rmin, outer_rmax;

    // 得到延拓后网格中一个四面体方位
    Point3d GetTetrahedronPosition(const Point3d &a, const Point3d &b, const Point3d &c, const Point3d &d);
};

class GeoGrid {
public:
    GRDECLData* data;
    // independent layer points
    LayerPoints* ind_layer_ps;
    GEO::Mesh mesh;

    GeoDelaunay delaunay;
    GeoDelaunay2d delaunay2d;
    // 网格尺寸
    double grid_size;
    // 网格生成策略：0-使用所有点直接划分，1-对于每层分别划分
    int strategy = 0;
    bool is_expand = false;

    LayerSurface* layer_surfaces;
    std::vector<int> tetrahedron_cells;
    std::vector<Point3d> mesh_points;
    std::vector<int> mesh_points_type;

    std::vector<int> tetrahedron_cells_expand;
    std::vector<Point3d> mesh_points_expand;
    std::vector<int> fault_color;

    ExpandInfo expand_info;

    GeoGrid(GRDECLData* data);

    ExpandInfo GetExpandInfo() {
        return expand_info;
    }

    double GetPillarZvalAtLayer(int layer, Line pillar);
    // 获取两点连线之间的中间网格点
    std::vector<Point3d> GetMiddleGridPoints(Point3d p1, Point3d p2);
    void setStrategy(int strategy) {
        this->strategy = strategy;
    }
    void SetGridSize(double grid_size) {
        this->grid_size = grid_size;
    }
    void SetExpandInfo(double expand_top, double expand_bottom, double expand_X_ratio, double expand_Y_ratio) {
        expand_info.expand_top = expand_top;
        expand_info.expand_bottom = expand_bottom;
        expand_info.expand_X_ratio = expand_X_ratio;
        expand_info.expand_Y_ratio = expand_Y_ratio;
    }
    // 计算和采样底面的边界和断层上的边
    void GetSurfaceBorderFault(std::vector<Point3d> &pts, std::vector<INDEX> &bdrs);
    // 根据底面信息三角化整个表面
    void SetSurfaceByBottomSurface(const std::vector<Point3d> &pts, const std::vector<INDEX> &bdrs);

    // 进行延拓
    void Expand();
};

#endif // GEO_GRID_H    